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The sign problem is a notorious problem, which occurs in Monte Carlo simulations of a system 
with a partition function whose integrand is not positive. One way to simulate such a system is to 
use the factorization method where one enforces sampUng in the part of the configuration space 
which gives important contribution to the partition function. This is accomplished by using con- 
straints on some observables chosen appropriately and minimizing the free energy associated with 
their joint distribution functions. These observables are maximally correlated with the complex 
phase. Observables not in this set essentially decouple from the phase and can be calculated with- 
out the sign problem in the corresponding "microcanonical" ensemble. These ideas are applied 
on a simple matrix model with very strong sign problem and the results are found to be consistent 
with analytic calculations using the Gaussian Expansion Method. 
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1. Introduction 



Monte Carlo simulation is an essential tool in the non perturbative study of a quantum field 
theory and has provided a wealth of information and insight in high energy physics. Its success lies 
in the fact that algorithmic importance sampling of a very small fraction of configurations gives a 
very accurate estimate of expectation values of physical observables. By stochastically constructing 
Markov chains of configurations, the system quickly thermalizes and one is able to sample the 
extremely small subset 3? of the configuration space that gives the important contribution to the 
partition function. ^ is determined by the competition of entropy with the Boltzmann factor of the 
action. For an observable ^ whose value on the configuration jx is ^[/x], the Markov chain sample 
{Mi jM2) • • • is constructed with samphng probability and (^)'s estimator is 

By choosing oc e^'^I'^'l , the above sum is easily performed numerically without terms whose size 
varies exponentially with the system size' . 

This approach breaks down in the case where the action of the system is complex and almost 
all configurations contribute a non positive term in the partition function. The sample is constructed 
by simulating a phase quenched model. One then mostly samples a subset of ^ whose over- 
lap with M is exponentially small with system size. The reason is that besides entropy and the 
Boltzmann factor of the real part of the action, the fluctuations of the imaginary part of the action 
become a determining factor in the suppression of configurations. Moreover, if one takes a simple 
reweighting approach to estimating (i^), one obtains 

where r[/i] is the imaginary part of the value of the action on /i. Then, besides the small overlap, 
one has to sum terms that oscillate wildly due to the exponential that appears in the numerator and 
the denominator. The oscillations are much stronger away from the stationary configurations of T 
which is usually what happens when we sample in Mq. 

An approach to studying such systems is the factorization method, originally proposed in 
used in [Q] and recently generalized in [||]. The key point is the selection of a maximal set of 
observables 

1^ = 1,... ,4 (1.3) 



which are strongly correlated with the phase e''^. The idea is that by solving saddle point equations 
for the minimum of the "free energy" associated with the distribution functions of these observ- 
ables, one can essentially determine the region The solutions are obtained by performing Monte 
Carlo simulations on selected subspaces of by constraining the values of these observables. Fur- 
ther calculations of any other observable 0' can be done by essentially sampling in ^ and their 



like e.g. in multihistogramming. 
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expectation values can be estimated without the phase factor appearing in (1.2). To be more spe- 
cific, if {xi ,X2, . . . ,x„} are the values^ of {^i , ■■ ,^n} that are the solutions to the free energy 
minimization equations, then one finds that 



^2,...^n (1-4) 

where (. . .)xi,x2,...,x„ are expectation values in a "microcanonical" system where are constrained 



to be equal to Xk- The right hand side of (1.4) has no complex action problem. The complex action 
problem has been reduced to computing the solution {xi, . . . ,x„} and this is greatly improved by 
factorizing the phase factor and taking advantage of its (hopefully) nice scaling properties. 

In this talk, the above general statements will be made concrete in a specific example. A 
pedestrian's approach is adopted for presenting the main ideas and the reader is referred to [||] 
for the details. A simple matrix model with very strong complex action problem is studied and is 
shown how to compute the set £ and the solution {xi , . . . ,x„}. Solving the equations that give the 
stationary points of T plays an important role in determining £. The nice scaling properties of the 
distribution functions are heavily used in the computation of the solution {xi, . . . ,x„}. 

Since the method can in principle be applied to any system with a complex action problem, we 
hope that our exposition will contribute to a successful application of the method to other interesting 
problems in lattice field theory and elsewhere. 

2. The Model 

Matrix models have been studied extensively in the context of non perturbative formulations of 
string theory and in the study of gauge/gravity duality. Monte Carlo simulations have contributed 
crucially in the study of the large N limit of supersymmetric matrix models [Q], in providing first 
principle evidence of gauge/gravity duality and in explaining the thermodynamics of certain black 
hole solutions in terms of microscopic string degrees of freedom 

The large limit of the IIB matrix model has been proposed as a non perturbative defini- 
tion of string theory. In this model, space-time emerges as the eigenvalue distribution of the bosonic 
matrices which makes possible to study the scenario of dynamical compactification of extra dimen- 
sions. This happens via spontaneous symmetry breaking of the rotational symmetry (SSB) of the 
model. Calculations using the Gaussian Expansion Method (GEM) support the realization of such 
a scenario [^. Monte Carlo simulations can play an important role in confirming those results 
from first principle calculations and in elucidating the mechanism that is responsible for SSB. The 
strong fluctuations of the phase e'^ favour length scales over which spacetime extends which are 
quite different than in the phase quenched model A simple matrix model that realizes the above 
scenario has been proposed in [||]. It is defined by the partition function 

Z = j dAdYdVQ-^''^''^ where Sb = ^ A^tr(A^)2 , Sf = -ij/^ (r^)«^ A^v^^ . (2.1) 
A^ in = 1,2,3,4) are N X N hermitian matrices, and iffa and i/A^ (a = 1,2, / = 1, . . . ,Nf) are 



A/^-dimensional row and column vectors. The actions (2.1) have an SO (4) symmetry, where the 



^More precisely are equal to the values of ffj^/ {ffk)o as it will be explained in the next paragraph. 
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bosonic variables transform as vectors and the fermionic variables transform as Weyl spinors. 
Integrating out the fermions, we obtain Z = JdA e^^'' Zf[A] , where Zf [A] = (det and ^ = F^A^ 
is a 2N X 2N matrix. The fermion determinant det & for a single flavor is complex in general but 
it turns out to be real for configurations with A4 = and that the phase of the determinant becomes 



stationary for configurations with A4 = A3 = []10|]. We take the \wge,-N limit with r = Nf/N fixed, 
which corresponds to the Veneziano limit. 

The order parameters of SSB of S0(4) are the expectation values of the eigenvalues of the 
"moment of inertia tensor" T^y = ^tr(A^Av) . These are ordered as Ai > A2 > A3 > A4 > and 
if their VEVs turn out to be unequal in the \wgQ-N limit, it signals the SSB of S0(4). GEM 
calculations show that S0(4) breaks down to S0(2) for all r > and Monte Carlo simulations 
are consistent with this result [^]. 

The "phase quenched model" is defined by 

Zo = j dA f'-^o W , Sq [A] = 5b [A] - A^f log | det &[A] \ . (2.2) 
and in this case SSB is absent. In fact (A„)o = 1 + 5 for all « = 1,2,3,4 where the VEVs ( • )o 



are taken with respect to (2.2). In the following, we will study the eigenvalues of T^y normalized 
to their phase quenched expectation values and we denote A„ = A„/ (A„)o. Deviation of from 
1 indicates a strong effect of the complex phase. In order to simulate ( ^ ) we rewrite it as Z = 
/ dA e^'^"''*! q'^W ^ where due to nonzero T[A] for generic configurations, the system turns out to have 
a very strong complex action problem. Due to the stationarity of the phase for lower dimensional 
configurations, the eigenvalues k„ are chosen for the application of the factorization method. 

We study the distribution function p(xi,X2,X3,X4) = (H^^i S{xk — Xk)) and the corresponding 
one in the phase quenched model p(*^^(xi ,X2,X3,X4) = {Y[t=i^{^k — ^k))o- By defining the "micro- 
canonical ensemble" Z^^ X2 X3X4= J dAe^^°Y[t=i ^{^k ^k) , we define the function w(xi,X2,X3,X4) = 
,.V2,JC3,JC4 and one finds that 

p(xi,X2,X3,X4) = ^p(°)(xi,X2,X3,X4)H'(xi,X2,X3,X4), (2.3) 

where C = (e''^)o. The minimum of ^(xi,X2,X3,X4) = — lnp(xi,X2,X3,X4) is the estimator for 
{(Ai), (A2), (A3), (A4)}. One has to solve the saddle point equations 

^/i°'(xi,X2,X3,X4) = -■^<I>(xi,X2,X3,X4), « = 1,2,3,4 (2.4) 



where f,^'' (xi , X2 , X3 , X4 ) = ^ In p (xi , X2 , X3 , X4 ) , <I>(xi , X2 , X3 , X4) = lim/v^oo p- In w (xi , X2 , X3 , X4) . 
The nice scaling properties of fli^^ ( ) and <5>(xi , X2 , X3 , X4 ) allow for extrapolations in N 

and x„. It is also important to note that the error in (A„) does not propagate exponentially in N. 



Looking for the full solution of ( |2.4| ) using Monte Carlo is a formidable task. Therefore, guided 
by the GEM results in we look for S0(3) and S0(2) symmetric solutions. For lack of space 
we discuss the S0(3) symmetric vacuum and the interested reader is referred to [^] for the S0(2) 
case. We take xi = X2 = X3 > 1 > X4 and define the reduced functions Pso(3) (■^'3') = P^°' {^T^,Xiy)7 
wso{3){x,y) = w{x,x,x,y), /so(3) -.^(x,^) = ^lnp^o'(3)(x,3;) and<I>so(3)(-^,3') =limw^oo ^lnwso(3)(-^>3')- 
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Figure 1: (Left) The function p-logwso(3)(-'Cj0.5) is plotted against x for N = S, 12, 16. The solid lines 
represent the function <I>so(3)(-x,0.5) obtained by extrapolation to = as described in [^. (Right) The 
function p-/so(3) x(^^^-^) plotted against x for = 16,32,64. The solid lines represent — ^<I>so(3)('^i0.5) 
obtained from the plot on the left. 
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S0(2) 
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factorization 


GEM 


factorization 


GEM 
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1.373(2) 


1.4 




1.151(2) 


1.17 


0.649(4) 


0.7 


(A4) 


0.59(2) 


0.5 


0.551(2) 


0.5 



Table 1: The results for the normalized eigenvalues (A„) for r = 1 obtained by the factorization method 
for the SO(3) and SO(2) symmetric vacua. The dash implies that the result should be the same as the one 
below in the same column due to the imposed symmetry. We also show the GEM results obtained at = 00 
inref. 



Equations ( 2A ) become 

^/sS(3),.(^'3') = -§^^so{3){x,y) , ^/^S(3),>.(^'3') = -§^^so(3){x,y) ■ (2.5) 

In the following we fix r = Nf/N = 1 and from the GEM results |^] we expect that the S0(3) 
symmetric solution is (x^^^ , y^^^ ) ^ ( 1 • 17 , 0.50) . We first fix j = y^,^^ =0.50 and solve for x. The 
results can be seen in fig. [T[ We find x = 1.152(3) which is consistent with x^^^ = 1.17. Next 
we fix X = X(,g^ = 1.17 and solve for y. We repeat the same procedure for the S0(2) symmetric 
vacuum and our results are summarized in Table [l[ 

Including more variables in the analysis is straightforward. If {X,X,Y,Z) is the absolute 
maximum of p{xi , X2 , X3 , X4 ), consider the "microcanonical ensemble" {•)x.x.¥,z and for an op- 

erator define p^(x) = ^^'^e^jf^;;^'''" . (^) = {8{x - ^))x,x,y,z where ^ = ^/(€?)o. Then 
Pei^) = r' Pfi\^)w&{x) where C = {e}^)xx.Y.z and w^(x) = ^^/^r^-;^'^^'''''^'''^ - Next we determme 

{&) from ^ff = -d.^ff{x), where = fjnpP{x) and ^ff{x) = limN^„^lnwff{x). If {0 
does not shift much from {^)x.x.y.z then the overlap problem for is not severe. Then one can 
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X X 

Figure 2: (Left) The function p-logw^ (x) is plotted against x for N = 6,8, 12. We also plot the function 
^(}{x) obtained by extrapolation to = oo as described in [||. The two solid lines represent the margin of 
error. (Right) The function -^j^\ogp'§\x) is plotted for = 8, 16,32. We also plot -^<I>^(x) obtained 
from the plot on the left. 



show that [||] 



X.X.Y.Z 



(2.6) 



{&'^)x.X,Y.^ 

i.e. correlation of and e''^ within {•)x,x,y,z is small. The advantage of this relation is that the phase 
factors out and one can calculate the expectation value of & within the microcanonical ensemble 
without the sign problem. 

We test the above relations by considering the observable = — ^L^^vtr[^M'^v]^- Asx<Cl 
the dominant configurations are approximately simultaneously diagonalizable [A^,Av] ~ and 
det ^ > 0. Therefore 0' can potentially have strong correlations with the phase factor. From fig. ^ 
(Left) we find that ^logw^'(x) approaches zero for ;c ^ as expected. From fig. ^ (Right) we 
find that the effect of the phase is to shift the estimate of {ff) by Ax = 0.07(3). On the other hand, 
the standard deviation of the distribution p^''(x) is estimated as a ~ 0.1 /N from the slope of the 
function plotted in fig. |2| (Right) around x ~ 0.92. This means that the deviation Ax is < 2 a for 
N < 16. Thus, the remaining overlap problem associated with this observable is practically small. 
This is consistent with the fact that we were able to reproduce the GEM result by constraining only 
the four observables X„ (n = 1,2,3,4). 



3. Conclusions 



In this work, we have applied the factorization method to the Monte Carlo study of a matrix 
model with strong complex action problem. It has been extended to include more than one observ- 
ables in order to eliminate the overlap problem and sample effectively the theory's configuration 
space. A maximal set £ of such observables has been determined that has significant correlations 
with the complex phase. By constraining sampling to configurations that lie in the neighbourhood 
of the solution to the minimization of the free energy of the distribution functions of the observables 
in r, all other observables can be computed without the complex action problem. The complex ac- 
tion problem has been reduced to solving the saddle point equations, which is difficult, but not 
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impossible if one takes advantage of the nice scaling properties of the factors in the distribution 
functions. Solving the equations for the stationary configurations of T played an important role in 
determining Z. 

Our results turn out to be consistent with analytical calculations using the GEM and confirm 
SSB of rotational invariance in the distribution of eigenvalues of bosonic matrices. This is related 
to the problem of dynamical compactification of extra dimensions in string theory which will be 
studied further in the ITB matrix model as well as in simpler matrix models with supersymmetry. 

The steps followed in this study are considered to be generic and we hope that a similar ap- 
proach can be apphed to many interesting systems. 
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